function [Xdraws_best] = func_likely_draws_current(IJitidx, IJchoiceset,IJcurrentFirst, Xdraws,Xmat,thetahat,y_logistic)
%This function takes the estimated coefficient estimates combined with the
%paramaeters of the random coefficient distribution and returns the 
%most likely draws for each individual based on evaluating the likelihood



% rename objects to match func_mixed_logit
sumidx=IJitidx(IJchoiceset==1,:);
sumidxC=IJitidx(IJcurrentFirst==1);
X=Xmat(IJchoiceset==1,:);
XC=Xmat(IJcurrentFirst==1,:);
draws=Xdraws(IJchoiceset==1,:,:);
drawsC=Xdraws(IJcurrentFirst==1,:,:);
theta=thetahat;
y=y_logistic;

beta = theta(1:size(X,2));
sigma = theta(size(X,2)+1:end);

sumdraws = zeros(size(draws,1),size(draws,3));
sumdrawsC = zeros(size(drawsC,1),size(drawsC,3));
for k=1:size(draws,2)
   sumdraws = sumdraws + sigma(k)*squeeze(draws(:,k,:));
   sumdrawsC = sumdrawsC + sigma(k)*squeeze(drawsC(:,k,:));
end

umat = (X*beta)*ones(1,size(draws,3)) + sumdraws;
expumatC = exp((XC*beta)*ones(1,size(drawsC,3)) + sumdrawsC);

%umatCfull = zeros(max(sumidxC),size(sumdrawsC,2));
umatCfull = zeros(max(sumidx),size(sumdrawsC,2));

umatCfull(sumidxC,:) = expumatC;

umatCfull = umatCfull(sumidx,:);

%XCfull = zeros(max(sumidxC),size(XC,2));
XCfull = zeros(max(sumidx),size(XC,2));

XCfull(sumidxC,:) = XC;

XCfull = XCfull(sumidx,:);

%drawsCfull = zeros(max(sumidxC),size(drawsC,2),size(drawsC,3));
drawsCfull = zeros(max(sumidx),size(drawsC,2),size(drawsC,3));

drawsCfull(sumidxC,:,:) = drawsC;

drawsCfull = drawsCfull(sumidx,:,:);


p1mat = max(min(exp(umat)./(1+umatCfull+exp(umat)),.9999999),.00000001);
p0mat = 1-p1mat;

pCmat = umatCfull./(1+umatCfull+exp(umat));

ymat = y*ones(1,size(draws,3));

probmat = p1mat.*ymat+p0mat.*(1-ymat);

Lmat = zeros(max(sumidx),size(draws,3));

for s=1:size(draws,3)
    Lmat(:,s) = accumarray(sumidx,probmat(:,s),[],@prod);
end

% pick out the most likely draw
[M,index_bestdraw]=max(Lmat,[],2);

% now figure out number of draws, and create an empty matrix to fill
Xdim=size(Xdraws);
Xdraws_best=zeros(length(probmat), Xdim(1,2));  % empty matrix for best draws
Xdim=Xdim(1,2); % number of random draws

% make the index_bestdraw long (essentially, copy it) 
indexlong=index_bestdraw(IJitidx,:);
indexlong=indexlong(IJchoiceset==1,:);
long=(1:length(indexlong))'; % create an index running from 1 to number of observations

% now fill the Xdraws_best matrix
for i=1:Xdim
    long1=i*ones(length(indexlong),1);
    bestindex=sub2ind(size(Xdraws),long,long1,indexlong);
    Xdraws_best(:,i)=Xdraws(bestindex);
end


end

